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ABSTRACT 

We analyze two time-dependent cluster cooling flow models in spherical symmetry. 
The first assumes that the intracluster gas resides in a static external potential, and 
includes the effects of optically thin radiative cooling and mass deposition. This cor- 
responds to previous steady-state cooling flow models calculated by White & Sarazin 

(1987) . Detailed agreement is found between steady-state models and time-dependent 
models at fixed times in the simulations. The mass accretion rate M is found either to 
increase or remain nearly constant once flows reach a steady state. The time rate of 
change of M is strongly sensitive to the value of the mass deposition parameter q, but 
only mildly sensitive to the ratio f3 of gravitational binding energy to gas temperature. 
We show that previous scaling arguments presented by Bertschinger (1988) and White 

(1988) are valid only for mature cooling flows with weak mass deposition [q < 1). The 
second set of models includes the effects of a secularly deepening cluster potential and 
secondary infall of gas from the Hubble flow. We find that such heating effects do not 
prevent the flows from reaching a steady state within an initial central cooling time. 

Subject headings: cooling flows — galaxies: clusters: general — hydrodynamics - 
intergalactic medium — X-rays: general 



Introduction 



X-ray and optical data strongly suggest that cooling accretion is taking place in more than 
half of galaxy clusters (for a review, see Fabian 1994). In non-cooling flow clusters, X-ray surface 
brightness profiles obtained from imaging data are usually well modeled by a single-temperature 
gas in hydrostatic equilibrium with a background isothermal potential (Jones & Forman 1984). In 
contrast, cooling flow clusters exhibit significant excess central emission compared to that derived 
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from a fit to an isothermal. This excess emission is thought to be due to gas losing its thermal energy 
and condensing out of the intracluster medium (ICM) at rates that often exceed 100 M yr _1 . Over 
a cluster lifetime, the central dominant galaxies in these clusters may accrete up to ~ 10 12 M Q of 
cooled gas. (We assume a Hubble constant of 50kms _1 Mpc _1 and that a cluster lifetime is 
comparable to a Hubble time.) 

X-ray surface brightness profiles of cooling flow clusters can be used to infer mass accretion 
rates M S urf which increase roughly linearly with radius from the centers of cooling flows (e.g., 
Stewart et al. 1984; Thomas, Fabian, Sz Nulsen 1987). If the flows are steady, this in turn implies 
that matter must be cooling and condensing out of the ICM over the full range of radii (> 100 kpc) 
for which cooling flow emission is observed. If there were no distributed mass deposition, then gas 
would be deposited only at the center. The resulting X-ray surface brightness profiles would be 
significantly more centrally peaked than those observed. 

X-ray spectra of cooling flows imply cooling rates M spcc within a factor of two of those inferred 
from the imaging data (e.g., Canizares, Markert, & Donahue 1988; Allen et al. 2000). This result 
is an independent confirmation of the accretion rates derived from imaging data. The dynamical 
accretion rate in a spherical system, M = 4irr 2 pu (where p is the gas density and u is the flow 
velocity) , is not observed directly. Although the dynamical and X-ray-derived accretion rates should 
be the same if the flow is in a steady state, the radial velocities in the accretion flow are typically 
on the order of tens of kilometers per second, well below both the velocity dispersion of a single 
galaxy and the spectral resolution of current X-ray instruments. 

Although the final repository for the cooling gas remains ambiguous, the evidence is quite 
strong that it is cooling and decoupling from the flow (Fabian 1994; White, Jones, & Forman 
1997; Peres et al. 1998; Markevitch et al. 1998; Allen et al. 2000; White 2000). Heating processes 
are unlikely to balance cooling since their functional dependences on density and temperature do 
not match those for radiative losses. Indeed, previous numerical investigations have shown that 
alternatives to large cooling accretion rates, such as heat conduction (Bregman & David 1988; 
Meiksin 1988), supernova heating, and drag heating by orbiting galaxies (Bregman & David 1989), 
are able to reproduce the data only for narrow ranges in the respective free parameters, if at all. 
Moreover, it is difficult to imagine a process other than cooling that can power an excess emission 
rate as high as 10 44 erg s" 1 without producing other noticeable signatures. Finally, most alternatives 
to cooling accretion fail to account for observed soft X-ray lines. 

Models that assume a steady accretion flow with mass deposition can reproduce the observed 
X-ray features of cluster cooling flows (e.g., White & Sarazin 1987, hereafter WS; Fabian 1994). 
WS calculated steady-state models with a simple mass deposition formation law, but such models 
provide no information on the evolution of cooling flows. 

In this paper we use a suite of time-dependent spherical models to assess the evolution of 
cooling flows in static and evolving gravitational potentials. In addition to providing predictions of 
the spatial structure and spectral properties of relatively relaxed cluster cooling flows, the spherical 
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model also has implications for their long term evolution. The time evolution of the accretion rate 
has consequences for the total amount of accreted material in cooling flows, and it can provide a test 
for the models when compared to X-ray observations of high-redshift clusters. Two free parameters 
likely to affect the time-variation of M in the spherical models are the mass deposition efficiency 
and the ratio (3 of gravitational binding energy to thermal energy in the gas. The parameter (3 
determines the shape of the gas distribution in clusters when the gas is in hydrostatic equilibrium. 

Self-similar models can offer insight into the evolutionary effects, and they provide a useful 
testing ground for fully time-dependent calculations (Chevalier 1987, 1988; Bertschinger 1989; 
Lufkin & Hawley 1993). However, they are limited in that they allow either a narrow range of 
initial conditions or a limited number of physical processes. As a result, such models have not 
been successful in attaining detailed agreement with the observations. Simple scaling arguments 
are potentially useful (Bertschinger 1988; White 1988), and we review these in § 2 below. Our 
numerical simulations of cooling flow evolution in static gravitational potentials are presented in 
§ 3. 

As a final point of inquiry, we examine the consequences of continued cluster evolution by 
investigating the effects of a secularly deepening cluster potential and continued accretion of gas 
from the Hubble flow. We present a simple physical argument in § 2.3 which shows that adiabatic 
compression of gas in a deepening gravitational potential does not inhibit cooling flows. Our 
numerical simulations of cooling flow evolution in an evolving gravitational potential are described 
in § 4 and compared to similar work by Meiksin (1990, hereafter M90). Our conclusions differ from 
those in M90, which found that cooling flows are slow to reach steady state and mass dropout rates 
were strongly reduced by time-dependent effects in an evolving gravitational potential. We do not 
find these effects in our simulations of cooling flows evolving in the same deepening gravitational 
potential as used in M90, which is consistent with the simple arguments we present in § 2.3. Our 
results are summarized in § 5. 



2. Physical Arguments 

Before proceeding to the numerical models, we first discuss some physical arguments regarding 
the temporal behavior of cooling flow clusters. These include definitions of the various physical 
regions in the flow, a review of analytic predictions of the time dependence of the cooling accretion 
rate M, and a discussion of whether adiabatic compression can extend the cooling time. These 
physical arguments form a basis for the discussion of the numerical models. 

The most important process in models for cooling flows is obviously the radiative cooling of 
the gas. Over essentially all of the region of the flow, this occurs on a timescale which is much 
longer than the dynamical time of the gas, and as a result the flows and compression which occur 
are relatively slow. We do not include thermal conduction or other diffusive processes in our 
calculations. The only other processes which affect the temperature of the gas are shocks and 
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adiabatic compression or expansion of the gas. In the absence of cluster mergers, the flow velocities 
are low and shocks are not important in our models. However, adiabatic compression of the gas 
will occur as a result of the cooling and inflow of the gas in the gravitational potential, or through 
the slow cosmological growth in the cluster potential. We discuss the effects of adiabatic heating 
on the cooling time scale of the gas in § 2.3 below. 



2.1. Timescales 

The general structure of cluster cooling flows in a static potential is determined by three 
characteristic timescales: the age of the cluster t agc , the cooling time t coo i and the dynamical time 
tdyn- The instantaneous isobaric cooling time is 

tco ° l ~ 2 pm p pA(T) ~ X [WkJ U * 10 24 ergscmS s" 1 g~ 2 J \W^F*) y "' ( j 

where p is the mean molecular weight in units of the proton mass m p (we assume p = 0.6 in this 
work), and n = pj ' pm p is the total particle number density. The function A(T) is the cooling rate 
such that p 2 A has the units ergs _1 cm -3 . We assume the cooling function of WS for half-solar 
abundances. The dynamical time is 

where r c is the cluster core radius and a c is the velocity dispersion of the cluster galaxies (also 
nearly equal to the initial sound speed in the gas prior to cooling). 

The cooling radius r coo i is defined as the point where the instantaneous isobaric cooling time 
equals the system age. Inside the cooling radius, the cooling time is less than the age of the system, 
and a cooling flow occurs. As long as the cooling time exceeds the dynamical time, the gas flows 
subsonically into the center, regulated by cooling. If the gas cools completely before reaching the 
center and the rate of mass drop out is not too large (q < 3), the flow generally passes through a 
sonic radius r s , where the flow speed equals the local sound speed (Sarazin & Graney 1991). For 
^cooi S> t S> r s the system is expected to be reasonably well-described as being in steady-state, 
an assumption which we test explicitly in § 3.3. The validity of the steady-state approximation at 
~ r coo i cannot be assessed except with time dependent models. 

The above definition for the cooling radius is straightforward, as well as familiar, but it may 
not be the most useful. Since the gas is flowing inward as it cools, we consider whether a more 
relevant timescale may be the integrated isobaric cooling time 

_ _5_ f O'dO' 
tint ~2Pj W)' 
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where 6 = kT/pm p . Typically, the integrated cooling time ti nt is about a factor of two shorter than 
the instantaneous cooling time i coo i- Consequently, gas may actually be taking part in the flow out 
to a considerably larger radius r- mt , where t- mt = t age . We comment more on this in § 3.2. 



2.2. M as a Function of Time 

The time-dependent calculations of §§ 3 and 4 solve explicitly for the accretion rate M as a 
function of time. However, it may also be possible to estimate M(t) based on the assumption that 
gas outside the cooling radius remains unaffected by cooling in the interior. Writing the dynamical 
accretion rate as M = ^i:r 2 pu, one might expect that the velocity can be approximated by the 
propagation rate of the cooling radius, so that M is interpreted as the rate at which gas is swept 
over by the cooling radius. This amounts to making the substitution u = dr coo \/dt\t cool . With 
this approximation, White (1988) showed that the cooling radius evolves with time as r coo \ cx 
where rj = [(1 — AfA)A r T — A r p]~ l : and where we have used the notation A x = (din /dlnx). The 
exponent of the time dependence in the accretion rate M is then given by 

A t M = (3 + A r p)rj -1 = 6 (4) 

If the gas is initially isothermal, then £ < (>) when A r p < (>) — 3/2. In other words, cooling 
flow accretion rates will increase (decrease) with time when the cooling radius is in a region where 
the gas density profile is shallower (steeper) than r~ . 

The gas density profiles of clusters of galaxies are often fit with the isothermal "beta" model. 
In this model, the gas density distribution is given by 



P(r) = Po 



1+1^ 



-3/3/2 

(5) 



where po is the central gas density (e.g., Jones & Forman 1984). We will sometimes give no = 
po/pm p , which is the corresponding total particle number density. If the galaxies in a cluster arc 
isothermal with a distribution given by the analytic King approximation to an isothermal sphere, 
then (3 is determined by the ratio of the velocity dispersion of the galaxies to that of the gas, 

(6) 

Typically, clusters have (3 2/3, and the density profile tends to steepen from A r p rj near the 
center to A r p r; —2 at large radii. Thus, there will be a transition radius r tran where £ = A t M 
changes sign from positive to negative. For cases where r coo i > r tra in on e would infer from equation 
(4) that M is decreasing. However, there appears to be some variation in the value of (3 from 
cluster to cluster, from 0.5 to about 1.2. This variation translates into a variation in the gas 
scale height relative to cluster core radii, and hence in the value of £ at the cooling radius. There 
is a correlation between cluster temperatures (both the gas temperatures at large radii and the 
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dynamical temperatures [velocity dispersions] of galaxies in clusters) and the overall slope of the 
X-ray surface brightness profiles (White 1991), implying that hot clusters (T gas >7x 10 7 K) tend 
to have density profiles scaling as p ~ r~ 2 at large radii, while cool clusters (T gas < 5 x 10 7 K) have 
shallower profiles. There is therefore a range of cases, with both increasing and decreasing mass 
accretion rates inferred via equation (4). In cases where imaging data imply a cooling flow, the 
slope of the X-ray surface brightness profile typically does not make a sudden inflection at either 
a core or cooling radius. Consequently, determining the transition radius from imaging alone is 
particularly difficult. Thus, in the absence of spatially resolved spectroscopy, there is some a priori 
ambiguity in the sign and magnitude of dM/dt. In § 3, we investigate this question for a plausible 



One obvious limitation to the analysis of § 2.2 is its failure to account for ongoing dynamical 
evolution of the cluster itself. As subclusters merge, and as matter continues to accrete from the 
Hubble flow, the cooling flow may be altered, or disrupted altogether. A 1-D calculation can take 
these processes into account only in the limit of quasi-static deepening of the cluster potential 
well. Observations (Bird 1994) and numerical simulations (Evrard 1990) of hierarchical structure 
formation show that the growth of galaxy clusters is only roughly approximated by a spherically 
symmetric deepening of a cluster's gravitational potential. If the changes in the cluster potential 
result from violent subcluster mergers, shocks may be driven into the gas. Such mergers cannot be 
easily modeled in a 1-D simulation. Nonetheless, following work by M90, we consider a simple model 
in which the growth of the cluster is described by a slow deepening of the cluster potential and the 
resulting adiabatic compression of the intracluster gas. Although this is a simplistic approximation, 
we can still address the question of whether such adiabatic compression is sufficient to change the 
qualitative nature of the steady flow compared to that of an isolated, static cluster, as has been 
suggested by M90. 

A simple physical argument shows that adiabatic compressive heating is not sufficient to bal- 
ance radiative cooling in intracluster gas. The point is that, whatever the cause of such compression, 
it is a reversible process and consequently does not change the entropy of the gas. Thus, the only 
changes in the entropy are those due to radiative cooling. Consider a parcel of gas initially at 
temperature Tj and density pi, with an adiabatic index 7. For the purpose of this argument, we 
approximate the cooling function as a power law, A oc T a . Adiabatic compression to density pf 
will heat the gas to a temperature of Tf = Ti(pf / ' pi) 1 ^ 1 . The cooling time scales as T 1_a /9 _1 , so, 
after compression, the cooling time is 



For 7 = 5/3 the effect of compressional heating is exactly balanced by the increased cooling rate for 
a critical value a = a* = —1/2. We note that a > —1/2 for almost all temperatures in the range 



range in (3. 



2.3. Adiabatic Compression 




(7) 
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10 6 <T< 10 9 , except for a narrow range between T ss 10 7 K and T«3x 10 7 K (e.g., Raymond, 
Cox, & Smith 1976; WS). Thus, adiabatic compression generally should decrease the cooling time 
for a given parcel of gas. This result is independent of the physical processes responsible for the 
compression. 



3. Time-Dependent Flow in a Static External Potential 

3.1. Assumptions and Equations 

The models in this section are designed to resemble those of WS for the purpose of assessing 
time-dependent effects in existing models with mass deposition. We therefore include radiative 
cooling and a mass deposition term, but ignore the effects of conduction and external heating, 
referring discussion of the latter to the published literature. Neglecting angular momentum, the 
corresponding spherically symmetric fluid equations are 

dp . P 9 2 

-dt + v^ ru) = - p (8) 

du ldP d$ _ 
dt p dr dr 

p|ln(Pp^) = -( 7 -l)p 2 A(T), (10) 

where P and T are the pressure and temperature respectively, p 2 A is the cooling rate in erg s _1 cm~ 3 , 
is the total gravitational potential and 

d_ = d_ d_ 

dt dt dr ' 

For the mass dropout rate p we adopt the cooling-time law of WS, 

P = qj^, (11) 

^cool 

where q is an efficiency factor of order unity and i coo i is defined in equation (1). This mass deposition 
law is based on the ansatz that thermal instabilities due to cooling lead to growing perturbations 
which leave the flow over a range of radii. (It is often expedient in numerical simulations to cut off 
the cooling function at some floor temperature Tfl oor , the value of which is arbitrary provided that 
the gas temperature does not reach Tfl oor unless the flow becomes supersonic.) 

We solve equations (8)-(10) using a time-explicit hydrodynamics code in spherical symmetry 
(Lufkin & Hawley 1993). The code is similar to existing second-order Eulerian codes except that 
it solves the fluid equations in Lagrangian mode with a globally conservative remap to the fixed 
grid. The code has been tested on a variety of problems, including self-similar cooling flows and 
cosmological accretion (for details, see Lufkin & Hawley 1993). 
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3.1.1. The Assumed Mass Deposition Efficiency q 

All of the cooling flow models considered here have a finite mass deposition efficiency coefficient 
q. If one assumes no mass deposition, the mass accretion rate is constant with radius in a steady- 
state flow, and the gas cools catastrophically at the center (WS; Meiksin 1988, 1990), i.e., the cooling 
is sufficiently strong for the gas to pass through a sonic point as it approaches the center. The 
catastrophe is prevented (i.e. the flow remains subsonic all the way to the center) for q > 3.4. The 
dynamical mass accretion rate M is also found to be roughly proportional to radius for q ~ 3 — 4. 
Although a smaller q cannot be ruled out observationally in all value as small as q = 

implies an X-ray surface brightness profile that is much too sharply peaked in the center. We 
therefore restrict our attention to the two cases q = 1 (below which the central density is too 
sharply peaked) and q = 4, and examine the effect of this variation on model flows. 

For these applications, the numerical grid covers the full dynamic range from scales of ~ 1 kpc 
near the center out to a radius of several Mpc, with a reflecting inner boundary (u — > as r — > 0). 
In the q = 1 case there is generally a sonic radius near r ~ 1 kpc in the steady-state models. 
The sonic radius is therefore unresolved by the grid spacing of ~ 1 kpc near the center in the 
time-dependent models. If an outflow boundary inner condition is used, the grid must resolve the 
region between the sonic radius and the origin, in order to have supersonic outflow at the inner 
grid boundary (Lufkin & Hawley 1993). Unfortunately, timestep constraints prevent the use of an 
outflow boundary in this case, owing to gas velocities in excess of 100 km s" 1 inward from the sonic 
radius. 



3.1.2. The Background Cluster 

The initial conditions for the time-dependent simulations are chosen such that the gas is 
isothermal and in hydrostatic equilibrium with the external potential <3?. Integrating the hydrostatic 
equation we obtain for the initial density profile 

p(r) t=0 = p exp{-[$(r) - $ ]/(fcWM™p)} , (12) 

where is the asymptotic gas temperature at large radii, and po an d 3>o are the density and 
potential at r = 0. The background potential for the cluster gas is assumed be that of a massive 
central galaxy $ 9 plus a smooth cluster component $ c . The cluster mass distribution is taken as 
a King approximation to an isothermal sphere: 

Pc (r) = p c0 [l + (r/r c ) 2 ]- 3/2 , (13) 
where r c is the cluster core radius. The potential for this distribution is given by 



In 



/r c + y/1 + (r/r c ) 2 



Mr) = -$ c o^ ; L ■ (14) 



r/r, 
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The central potential $ c o = 4irGp C Qr^ and is related to the line-of-sight velocity dispersion a c 
by $c0 = 9cr 2 . We assume a c = 1054 km s -1 (<E> c o = 10 17 erg g _1 ) for the hot cluster runs, and 
a c = 527 km s -1 for the cool cluster runs. The corresponding cluster masses inside of 250 kpc arc 
10 14 and 5 x 10 13 M Q , respectively. 

We note that our adopted potential has a finite depth, whereas a true isothermal is infinitely 
deep. Consequently, the density profile for the gas becomes sufficiently shallow at large radii that 
the total X-ray luminosity diverges. This actually follows from assuming that the gas at large radii 
is hydrostatic and isothermal, which cannot be true on spatial scales where the sound crossing time 
is greater than a Hubble time (r ~ lOMpc). In the evolving cluster models of § 4, the outer parts 
of the cluster are in free fall from the cold Hubble flow, with an accretion shock at a radius of 3-10 
Mpc. Because the cooling flow is confined to the inner 100 kpc or so, the solution is not sensitive 
to the structure of the cluster much beyond this radius. We therefore quote total luminosities 
within one half and two Mpc. Furthermore, because gas densities far from the center are not well 
constrained observationally, we omit discussion of the structure of the ICM at very large radii. 



3.1.3. The Central Galaxy 

All of the models here and below assume the presence of a massive, stationary central galaxy. 
However, in choosing a form for the galactic potential we are guided by the desire to extend these 
models to higher dimensions, with the possibility of allowing a galaxy to orbit in the cluster (e.g., 
Lufkin, Balbus, & Hawley 1995). The primary considerations are computational efficiency and 
accuracy when the galaxy's position changes with time. The potential should be smooth and fully 
resolved, with — > as r — > to avoid grid noise. It should also be expressible in terms 
of simple analytic functions. A potential of the form of equation (14) would be one possibility. 
However, the mass in this model diverges logarithmically, and it is generally preferable to have a 
galaxy potential corresponding to a finite mass. While it is tempting to truncate an analytic King 
model at some finite radius, this can generate anomalous sound waves for an orbiting galaxy on a 
finite difference grid. 

A simple form which satisfies our numerical criteria is Plummcr's model (Binney & Tremaine 
1987): 

$ 9 = cD 90 (l+r 2 /r|)- 1/2 , (15) 

where r g is a characteristic radius. As it stands, the density distribution which gives rise to equation 
(15) does not approximate real galaxies well. However, we have found that a superposition of 
suitably weighted components can give a good approximation to a King model. For example, the 
addition of two potentials, one for the center and one for an extended distribution, can be written 

as 



0--v) + v 



(l +r 2/ r 2)l/2 (1 +r 2/ r 2)l/2 



(16) 
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where < y < 1 and r t is the characteristic radius of the extended component. We have chosen to 
specify the potential directly, rather than starting from the underlying density distribution. This is 
because the gas interacts with the galaxy only through gradients in the potential. The analytic form 
we have chosen contains functions which can be calculated at high speed, whereas transcendental 
functions (such as power laws, logarithms and arctangents) can increase the run time significantly. 
The square root function is an exception. It uses a divide-and-average algorithm, which converges 
very fast. A tabulated function, although not particularly expensive, can be rather unwieldy, 
especially when a galaxy's position changes with time. 

At large radii, (r » r t ), the potential tends to that of a point mass, so that the total mass has 
the finite value 

M g = Klz%+jgdggg . (17) 

Gr 

For y = 0.15 and r\ = 40r g , the gradient of this potential closely approximates that of an analytic 
King model truncated at r = 20r g . 

We note that in equations (13) and (14), r c is a true core radius (i.e., it corresponds to the 
radius at which the surface density falls to one half of its central value) . In the model galaxy chosen 
here, this occurs at a radius of about 0.63r 9 . Figure 1 shows g = —dQ g /dr for a King model in 
which the central potential is <I> 9 o = 7<r^ (solid line; see Binney Sz Tremaine 1987), compared with 
truncated (short dash) and nontruncated (dotted) King approximations to an isothermal sphere. 
These are plotted along with a model galaxy with a potential of the form of equation (16), for 
y = 0.15, r g = 1.58r c , r t = 10r c (long dash). The vertical scale is in arbitrary units, and can be set 
either by the central potential or the total mass of the galaxy M g . 

In the simulations below, the core radii r c are 250 kpc for the cluster and 4.41 kpc for the 
central galaxy (r g = 7 kpc). Assuming a total mass for the galaxy of 2 x 10 12 M Q implies a line-of- 
sight velocity dispersion of 288 km s _1 at the center and 251 km s" 1 averaged over a circular aperture 
subtending a radius of 10 kpc, in the manner of Bailey &: MacDonald (1981). The composite density 
distribution can also be reasonably approximated by an NFW profile (Navarro, Frenk, & White 
1997) with a scaling radius of r s ~ 200 kpc. 



3.2. Results 

Table 1 lists the input parameters and some of the calculated global characteristics for ten 
simulations. Each of five model clusters are evolved for 15 Gyr for q = 1 and q = 4. The first 
model, (runs 1 and 2) is taken to be typical of a strong cooling flow in a hot cluster, and is discussed 
in detail in §§ 3.2.1. The second model (runs 3 and 4) corresponds to a relatively cool cluster and 
is discussed in § 3.2.2. The remaining three cluster models, also discussed in § 3.2.2, are used 
to explore a range of values in j3. The results of runs 1-4 are plotted in Figures 2-9, and are 
compared to steady-state calculations in § 3.2.3. 
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The physical assumptions, as stated above, are identical to those in WS, with the exception of 
the form of the external potential. Self-gravity is neglected for the static models, but is included 
in the evolving cluster models in § 4 below. Because the cooling flow is primarily pressure-driven, 
the inclusion of self-gravity has little effect. 

3.2.1. Hot Cluster Runs 

Shown in Figure 2 are the density, temperature, mass accretion rate and cooling times as 
functions of radius at various times in the evolution, for the hot cluster run with q = 1. The initial 
conditions, given in Table 1, are chosen so that the models have reasonable global properties at 
t = 10 10 yr in the evolution. In discussing the properties of the models, we use an age of 10 10 
yr as a fiducial time for comparison to observed clusters. The solid curves in the figures refer to 
this fiducial age. More generally, the dotted, short-dashed, long-dashed, solid and dot-dashed lines 
correspond to ages of 0, 3, 6, 10 and 15 Gyr, respectively. The curve at 15 Gyr (beyond the fiducial 
age) is included to show that a steady state has been reached. 

Fi gure 3 shows the time evolution of both the instantaneous and integrated cooling radii, r coo i 
and r- mt . We also give the dynamical mass accretion rate M = —A.nr 2 pu as determined at each 
of these cooling radii. The X-ray-derived mass accretion rate (or cooling rate) Mx is determined 
by assuming that all of the X-ray luminosity from within the cooling radius is due to gas cooling 
within that radius. This is the standard assumption made in analyses of the X-ray surface brightness 
profiles to derive the total cooling or accretion rates M sur f (e.g., Thomas et al. 1987). (Recall that 
Mx = M only for steady state flows.) As expected, there is close agreement between M and Mx 
at each radius, although the value of M at r- mt is nearly twice that at r coo i. The panel at lower 
right shows the value of £ = AtM as derived from local gradients in the simulation (eq. 4). We 
find £ ~ 0.29 at r int and £ ss 0.59 at r coo i. The value at r coo i compares reasonably to a value of 
AfM ~ 0.34 measured directly from the simulation by least-squares fitting of a straight line to 
log M(i) as measured at r coo i and dumped every 10 8 years. 

We determine the transition radius (where £ = in eq. [4]) by examining the numerical 
results directly: r tra n = 234 kpc at t = 10 10 yr. The measured cooling radius at this time is 
^cooi = 100 kpc < r tran , so we would infer from equation (4) that the accretion rate is increasing, 
as Figure 3 confirms. There are several sources of error in the estimate of r coo i. The grid spacing 
is about 5 kpc in this region, although this can be reduced if necessary by running the problem at 
higher resolution. More serious is the assumption that the relevant characteristic radius is r coo i, 
not ri n t ( see ec i s - 1 an d 3). At t = 10 10 yr in the q = 1 case, the location where t; n t = t a ge is 
Hnt = 151 kpc < r tra n- Thus we see that the scaling relation of § 2.2 is confirmed in this case for 
both measures of the cooling radius. 

Several other features are noticeable from the figures. Gas in the core has cooled completely 
between 6 and 10 Gyr, with little change in the flow thereafter. The initial cooling time at the center 
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is 10.8 Gyr, but the integrated isobaric cooling time is 6.43 Gyr. This time period corresponds to 
the phase in which the cooling radius and mass accretion rates slow to a steady rate of increase 
(Figure 3). The particular shape of the temperature profile is determined by the cooling function 
and the shape of the potential, while the density and accretion rate profiles are nearly featureless. 
Because the sonic radius is unresolved, the innermost grid zones reflect low-amplitude sound waves 
into the flow. These are visible in the accretion rate at 15 Gyr, but do not significantly affect the 
flow. 

One might expect an inflection point in the density profile as the center begins to cool, but 
this is not observed. The entire inner region cools almost simultaneously, but isobarically, resulting 
in a smooth transition to the steeper density profile seen at later times. 

The hot cluster q = 4 solution (run 2) is shown in Figures 4 and 5. This model has a higher 
rate of mass drop out at larger radii. In this case, the initial density required to recover an accretion 
rate of about 300 M yr _1 at t = 10 Gyr is sufficiently high (no = 0.033 cm" 3 ) that the center cools 
in about 5 Gyr. After this time we find that M is slowly decreasing, even though both cooling 
radii (r coo i = 74kpc; r- mt = 128 kpc) lie well inside r tra n = 354 kpc. Moreover, the logarithmic 
rates of change of Mx measured at both ri nt and r coo i are approximately equal £ « —0.45, whereas 
the value predicted by equation (4) goes no lower than £ = —0.2 for this time frame. The value 
of T) = d\nr coo \/dlnt is also discordant; we measure r] ~ 0.5 from the simulations, whereas the 
density and temperature profiles would suggest a value of r/ ~ 1.0. Thus, it appears that the 
scaling argument detailed in § 2.2, which assumes q = 0, breaks down for large values of q. We will 
comment more on this below. Because there is no sonic point in this model, there is no need for 
high resolution in the innermost regions of the flow. 



3.2.2. Cool Cluster Runs 

We perform a second pair of simulations in which all the input parameters are the same as 
those above except that the cluster has a central velocity dispersion of a c = 527 km s -1 and the 
asymptotic gas temperature is = 3 x 10 7 K. These parameters imply (3 = 0.6 (run 3 and 4 in 
Table 1), compared with f3 = 0.75 in the hot cluster simulations. The results for runs with q = 1 
and q = 4 are shown in Figures 6-9. The same general comments regarding the hot cluster models 
also apply in this case. The q = 1 case shows an accretion rate M that is increasing with time, 
in agreement with the arguments in § 2.2, while the q = 4 case shows a nearly constant M even 
though the predicted £ is clearly positive near the cooling radius and nowhere negative. 

There is some variation in the particular values of £ and rj compared with the hot cluster runs. 
To see if these quantities correlate with the value of [3 in these models, we ran a series of cool cluster 
models in which only the initial temperature was allowed to vary. All other external parameters 
are the same as those in runs 3 and 4. These are listed as runs 5-10 in Table 1, and the derived 
values of i] and £ (taken at r coo i at t = 10 10 yr) are plotted versus f3 in Figure 10. Other than the 
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extreme case of (3 = 0.36, (which is included only to show that we recover the expected limiting 
case £, rj — > oo as (3 — > 0), the variation with (3 is slight compared to the dependence of £ and r] on 
q. In every case the values of £ and r\ derived by measuring the time derivatives directly from the 
models are much closer to the expected values when q = 1 than when q = 4. We must therefore 
conclude that the scaling argument of § 2.2 applies only in the limit of small q. 

The reason for the inability of equation (4) to predict the flow behavior for large q is not 
immediately clear, although there are several ways in which the flow may be affected by a large 
value of q. One is that the rapid mass deposition can significantly reduce the gas density in the 
interior, thereby reducing the pressure, so that gas slumps to the center from radii well beyond 
the cooling radius. Unfortunately, there does not appear to be a well-defined position at which to 
evaluate £ from observed density and temperature profiles. This is because q is not known in real 
clusters, and it may be that the relevant cooling radius should be redefined again to account for the 
high rate of mass deposition in this model. Another possible cause for concern is that if the mass 
deposition formula (eq. 11) is applied throughout the flow, it implies that matter can condense out 
even if the local cooling time exceeds the age. To guard against this, we ran models 5-10 with a 
cutoff to the mass deposition (q — ► 0) for t coo i > t agc . We also repeated the first four runs with this 
modification, but there were no differences in the results other than a leveling-off in the dynamical 
accretion rate at the cooling radius. 



3.3. Correspondence to Steady-State Models 

We have seen that the dynamical mass accretion rate is a good approximation to the accretion 
rate one would infer, under the assumption of steady flow, from X-ray observations of the model 
clusters described above. A further test of the assumption of steady flow is to compute steady-state 
models corresponding to a fixed time in the evolutionary models. The equations for steady flow are 
given by WS, and are obtained by setting all partial time derivatives to zero in equations (8)-(10). 
The resulting ordinary differential equations are then solved numerically subject to appropriate 
boundary conditions. Because the time-dependent calculation used here and the ODE solver arc 
independent, this is also a further test of the hydrodynamics code. 

There are a number of possible ways to specify the three boundary conditions on the steady- 
state models. Following WS, we choose to fix the age (t agc = 10 10 yr), and values for the temperature 
T| rcool and the dynamical mass accretion rate M\ Tcool at the cooling radius. For direct comparison 
with the time-dependent models, we can take the values from Table 1. The solution to the steady 
flow equations is plotted in Figures 11 (hot cluster) and 12 (cool cluster). These curves can be 
compared directly to Figures 2 and 4 at t = lOGyr (solid lines). The agreement is very close, with 
only a small departure near r = in the models with q = 1. This difference is not unexpected 
given the artificial boundary conditions u | r= o = imposed in the time-dependent calculations on 
flows with sonic radii. 
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4. Time-Dependent Flow in a Static Evolving Potential 

One possible limitation to the preceding models is the neglect of cluster evolution. This can 
consist of dynamical evolution in the form of continued mergers or two-body encounters which 
act to deepen the cluster potential, or of continued infall of gas from the Hubble flow. The gas 
itself may contribute to the change in the total potential. M90 included all of these effects in a 
series of simulations in spherical symmetry. Two classes of models were considered: one with a 
centrally dominant galaxy (because all cooling flows have one) and one without a central galaxy. 
We reexamine the effects of cluster evolution on cooling flows and use the same cluster models as 
M90 so we may compare results. We are interested in whether the flows reach a steady state with 
or without mass deposition. 

4.1. Model Characteristics 

In the models of M90, the cluster potential is as described above in equation (14), with an initial 
velocity dispersion of about 500 km s -1 . The depth of the potential is assumed to increase gradually 
with time until, at the end of each run, the velocity dispersion is about 1000 km s -1 . This is done 
to mock the dynamical evolution of the cluster from a redshift of z = 1.5 to z = 0, corresponding 
to a time span running from 3.3 Gyr to 13.1 Gyr (for f2o = 1 and Hq = 50kms _1 Mpc" 1 ). The 
gas was assumed to trace the dark matter initially and the central gas temperature put as much 
energy per unit mass in the gas as in the dark matter. The initial temperature profile was assumed 
to be adiabatic. 

Models without a central galaxy assumed that the initial gas velocity follows the same profile 
as that for collisionless infall onto a region of overdensity equivalent to that given by the total 
cluster mass interior to each point. The gas is therefore at rest only at r = and at the turnaround 
radius far outside the cluster, with infall interior to the turnaround radius. Initially the gas in the 
cluster is thus not in hydrostatic equilibrium, but instead falls inward, then bounces adiabatically. 
The gas in the center relaxes on a dynamical timescale, which is much less than the cooling time, 
so the precise initial conditions do not seriously affect the subsequent evolution of the flow. M90 
found that gas will cool catastrophically at the center of such flows if q = 0, but that for q = 1/2 
the flow is stabilized against runaway cooling. (Note that M90 employs a different definition of q, 
and the value of q = 1/2 in M90 corresponds to q = 5/6 in our notation.) 

We are primarily interested in the models of M90 which contain a central galaxy in the initial 
conditions. In this case, the initial velocity profile is the same as in the case without a central 
galaxy except that u = for r < 250 kpc, resulting in a discontinuous initial velocity profile. 



-15- 



4.2. Numerical Simulations 

Figure 13 shows our simulations using the PLPC code (Lufkin & Hawley 1993) for the fiducial 
model of M90 with a central galaxy and without mass dropout (q = 0). All of the parameters given 
in M90, including the potential for the central galaxy and the self gravity of the gas, are duplicated 
here. (We use a different cooling function, which we comment on below.) Up until the second time 
frame, our results and M90 agree in detail. We note that the discontinuous velocity profile results 
in an impulsive squeeze on the initially hydrostatic atmosphere interior to r c . This strong pressure 
wave travels into the center, bounces off the reflecting inner boundary, and travels back out into the 
infalling gas beyond r c . This is evident in both the velocity and temperature profiles at t = 4Gyr 
(dashed line in Figure 13). We find, using fine time resolution in the data dumps, that two strong 
bounces occur before the gas relaxes at the center. (The relevant plots in M90 also show evidence 
for these bounces.) The run begins at t = 3.3 Gyr, when the central instantaneous isobaric cooling 
time is 9.0 Gyr. 

For time slices later than 4.0 Gyr, our results differ from those in M90. We find that the 
gas cools in less than an initial cooling time, with an ensuing cooling catastrophe at the center 
(t = 9.9 Gyr, dot-dash line of Figure 13). For the same model, M90 found no cooling catastrophe. 

We have made a number of tests to search of the source of this discrepancy with the results 
of M90. First, the calculations were repeated using exactly the same grid of 150 zones as used in 
M90: 

n = 1.88kpc(i-l) i < 21 

n = 1.046r;_i % > 21 . 

For comparison, we also made plots which were 1-2-1 smoothed (as in M90), but detected no 
variation from those shown in Figure 13. Second, we considered the effects of differences in the 
cooling function in the two sets of models. The cooling function adopted in M90 is that of Gaetz 
& Salpeter (1983), whereas we have used that of WS, which is somewhat stronger. We therefore 
ran the case with a central galaxy with the strength of the cooling function reduced by 20% to 
compensate. The result was virtually unchanged. This is not surprising, since the cooling rate is 
proportional to the square of the density. It therefore appears unlikely that the disparate results 
are caused by small differences in the cooling functions. Finally, we have run the fiducial model 
with two other hydro codes, the direct-Eulerian MONO scheme of Hawley, Smarr, & Wilson (1984) 
and a modified version of the implicit scheme of Ruppell & Cloutman (1975). Although the exact 
time of cooling collapse varies by about 20% (this measure is very sensitive to the value of the 
central density, which in turn is sensitive to the numerical diffusion as the impulsive compression 
waves bounce off the center early in the simulation), the density, velocity and temperature profiles 
agree in detail with the calculation of Figure 13. 

One possible reason for this discrepancy is the grid noise (zone-to-zone variations in the fluid 
variables) present in the numerical calculations of M90. This could cause sound waves to be 
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amplified or converted into heat energy, perhaps via artificial viscosity. We also note that the 
crests in the temperature and velocity profiles are sharper in our calculation than in the figures 
of M90. This is true even when our model output is smoothed identically to that present in 
M90. This suggests that the features in the plots in M90 are spread over several grid zones, which 
might indicate that a significant level of numerical diffusion was present. Such diffusion may have 
transported heat energy into the core. 

We find that a cooling catastrophe occurs either with or without the central galaxy. We 
do reproduce the result in M90 that the catastrophe occurs sooner in the case without a central 
galaxy even though the initial cooling time is longer. As suggested in M90, this is likely due to the 
homologous collapse of the constant-density core. 

When mass deposition is included, the flow reaches a steady state on an initial cooling 
timescale. The details in the flow variables appear to be affected by transient features left over 
from the discontinuous initial conditions (q = 1 and 4; Figures 14-17). In particular, there are dis- 
continuities in the position of the cooling radius due to the shock waves that bounced off the center 
at the beginning of the runs. As expected, the temperature beyond 100 kpc increases gradually 
due to continued shock heating from secondary infall. 

The agreement between the cooling and dynamical accretion rates is also quite close, as in the 
static models of § 3. Thus, we conclude that flows with mass deposition do reach a steady state 
when cluster evolution and secondary infall are included in the spherical model. For both runs, M 
stays constant to within a factor of two after the flows reach a steady state. Evolutionary effects 
in spherical symmetry therefore do not appear to affect the time evolution of the mass accretion 
rate relative to that found for static clusters. 

5. Discussion and Conclusions 

We have solved for the time-dependent behavior of cooling gas in a variety of spherically 
symmetric cluster models, with particular emphasis on the evolution of the mass accretion rate. 

We find that the steady-state approximation is valid after the initial onset of cooling at the 
center for cluster flows in both static and evolving external potentials. This result is insensitive 
to the inclusion of self-gravity. For the models considered here the accretion rate either increases 
or stays about the same with time. While the rates of increase or decrease that we see in the 
simulations would be difficult to infer from imaging observations, spatially resolved spectra can be 
a useful diagnostic (e.g., Wise & Sarazin 1993). The difference arises from the fact that while the 
age is an assumed parameter in the models, it can be measured by calculating the cooling time 
at the radius of extent of soft X-ray lines. Moreover, the cooling rates derived from imaging are 
more model-dependent, than those derived from spectra, which can in principle be obtained from 
a measure of the mass cooling through a single line. 
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A decreasing accretion rate with time is seen only for cases where the temperature of the gas 
corresponds closely to the cluster velocity dispersion ((3=1). In a number of clusters, especially 
poorer ones, the gas may be hotter by a factor of two than the virial temperature measured from 
galaxy motions. Thus, although we cannot rule out a decreasing mass accretion rate on the basis 
of these experiments, nearly constant or increasing accretion rates are favored for most clusters, 
unless they have been strongly affected by mergers. In the absence of mergers, the spherical model 
with mass deposition would predict that cooling flows were not much more vigorous in the recent 
past than they are today. Attempts to include additional evolutionary effects via secular deepening 
of the cluster potential well and continued accretion from the Hubble flow in spherical symmetry 
do not alter these results. However, it is likely that cluster mergers can have a strong effect on 
cooling flows (e.g., McGlynn & Fabian 1984; Gomez et al. 2000). Further study of the effect of 
cluster evolution on cooling flows will require numerical simulations in three dimensions with high 
spatial resolution and including cooling. 
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Fig. 1. — Plot of gravitational acceleration g versus radius for four model galaxies. The dotted 
line is a King approximation to an isothermal sphere (p ~ r _3//2 at large r). The short-dashed 
line corresponds to the same distribution truncated at 20 core radii. The long-dashed line gives 
the acceleration due to the superposition of two Plummer's models (derived from equation 16), 
with the parameters as stated in the text. The solid line is for a true King model with the central 
potential equal to la 2 c . 

Fig. 2. — The fully time-dependent solution for the hot cluster q = 1 model (run 1 in Table 
1). The four panels show total number density, temperature, dynamical mass accretion rate, and 
instantaneous isobaric cooling time as functions of radius shown at t = (dotted line), 3 (short 
dash), 6 (long dash), 10 (solid) and 15Gyr (dot-dash). 

Fig. 3. — Comparison of accretion rates for the q = 1 model of Figure 2. The top left panel 
shows the time evolution of the cooling radius determined from the instantaneous isobaric cooling 
time (eq. 1, solid line) and the integrated isobaric cooling time (eq. 3; dotted line). The upper 
right and lower left panels show the dynamical mass accretion rate M = —A7rr 2 pu and the cooling 
mass accretion rate derived from the X-ray emission Mx, respectively. The solid (dotted) curve 
corresponds to taking the accretion rate at the radius where the instantaneous (integrated) cooling 
time equals the age. The panel at lower right gives the estimated logarithmic time derivative of 
the mass accretion rate as a function of cooling radius (eq. 4) at the end of the run. 

Fig. 4. — The time-dependent solution for the hot cluster q = 4 model (run 2 in Table 1). The 
notation is the same as Figure 2. 

Fig. 5. — Same as Figure 3, but for the hot cluster q = 4 model (run 2 in Table 1). 

Fig. 6. — Cool cluster model with q = 1 (run 3 in Table 1). The notation is the same as in Figure 
2. 

Fig. 7. — Cool cluster model with q = 1 (run 3 in Table 1). The notation is the same as in Figure 
3. 

Fig. 8. — Cool cluster model with q = 4 (run 4 in Table 1). The notation is the same as in Figure 
2. 

Fig. 9. — Cool cluster model with q = 4 (run 4 in Table 1). The notation is the same as in Figure 
3. 

Fig. 10. — Top panel: the logarithmic rate of change of the cooling radius with time as a function 
of (3 for flows in a static potential. Bottom panel: the logarithmic rate of change of the dynamical 
mass accretion rate with time as a function of (3. Filled figures correspond to models with q = 4 
and open figures to q = 1; The triangles represent values that one would estimate by applying 
equation (4) to ideal observations of the model clusters, while the squares are the values obtained 
by direct measurement from the simulations. Here, the failure of equation (4) results from two 
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separate causes. In the q = 1 models the cooling flows are too young for complete cooling in the 
center, and in the q = 4 models the strong mass deposition results in structural changes that reach 
far beyond the fiducial cooling radius. Because the value of q in real clusters is not known, we must 
conclude that it is very difficult to measure reliably the rate of change of M directly from X-ray 
images. 

Fig. 11. — Steady-state cooling flow solutions for the hot cluster runs in Figures 2 and 4 obtained 
using the method of White & Sarazin (1987). The four panels show total number density, temper- 
ature, mass accretion rate and instantaneous isobaric cooling time as functions of radius for q = 1 
and q = 4. 

Fig. 12. — Same as Figure 11 except for the cool cluster runs of Figures 6 and 8. 

Fig. 13. — Our simulation of the fiducial run of M90 with a central galaxy and q = 0. The 
evolution times are 3.3 Gyr (initial conditions; dotted line), 4.0 (short dash), 7.0 (long dash) and 
9.9 Gyr (dot-dash). A cold core has developed in the innermost grid zone by 9.9 Gyr, preventing 
continued numerical evolution. 

Fig. 14. — Same as Figure 13 but for q = 1. The dot-dashed and solid lines correspond here to 
t = 10.0 and 13.0 Gyr. 

Fig. 15. — The evolution of cooling radius parameters in our simulation of the the q = 1 M90 model 
(same notation as Figure 3). The structure in the time-dependence of the mass accretion rate is 
due to transient sound waves and shocks resulting from non-static initial conditions. Despite these 
effects, the dynamical and cooling accretion rates trace each other quite well, indicating that the 
flow has reached a steady state. 

Fig. 16. — Same as Figure 14 but for q = 4. 

Fig. 17. — Same as Figure 16 but for q = 4. 
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